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The spUtting of e'^'-'^^^^ into a single product of e'''* and e'^^ results in symplectic integrators 
when A and B are classical Lie operators. However, at high orders, a single product splitting, with 
exponentially growing number of operators, is very difficult to derive. This work shows that, if 
the splitting is generalized to a sum of products, then a simple choice of the basis product reduces 
the problem to that of extrapolation, with analytically known coefficients and only quadratically 
growing number of operators. When a multi-product splitting is applied to classical Hamiltonian 
systems, the resulting algorithm is no longer symplectic but is of the Runge-Kutta-Nystrom (RKN) 
type. Multi-product splitting, in conjunction with a special force-reduction process, explains why 
at orders p = 4 and 6, RKN integrators only need p — 1 force evaluations. 



I. INTRODUCTION 

The approximation of c'^(^+^) to any order in h via a single product decomposition 

i 

is the basis of the splitting method^'^'^ for solving diverse classical^'^'^'^'^ i^'^^i^^'^^ , quantum mechanicali^ii^ii^ii^, 
and stochastic evolution equations^^'^^. The resulting algorithms are then symplectic, unitary and norm-preserving, 
respectively. However, at high orders, a single product splitting of the form ([!]) is very difficult to derive and the 
number of operators grows exponentially. This is easy to see: the error terms of the decomposition are governed by 
the fundamental Baker-Campbell-Hausdorff formula, 

^fiA^hB ^ ^h(A+B) + ^h^[A,B] + ^h='aA,[A,B]]-[B,[A,B]]) + --- ^2) 

where -B] = AB — BA is the commutator bracket. Starting with the first error term [A, successive higher 
order error terms are generated by further commutators of A and B, such as [A, [A,B]], [B, [A,B]], [A, [A, [A,B\\\ 
and [B, [A, [A, i?]]], etc.. This rapid doubling is moderated somewhat by the Jacobi identity and the special form 
of operators A and -B, but the growth of these error terms is essentially exponential. This implies that high-order 
symplectic integrators must require many force-evaluations to eliminate these many error terms. For symplectic 
integrators, the minimum number of force-evaluations necessary at orders 4, 6, 8, 10 has been found empirically to be 
3 (Ref-ii^i^iia), 7 (RefJ^), 15 (Ref^) and 31 (Ref.^O), respectively The doubling trend is clearly visible and seems to 
be precisely twice plus one. 

While symplectic integrators have many desirable conserving properties, such as the exact conservation of Poincare 
invariants, they are not immune from the irreversible phase error o^^'^^i'^^'^'^ that directly affects the accuracy of 
trajectories. Since the accuracy of the trajectory is of paramount importance, regardless of how well some Porincare 
invariants are preserved, it is less clear then, given the same number of force-evaluations, that a 2n-order symplectic 
integrator should be preferred over a 2n -I- 4 or 2ri -I- 6 order non-symplectic integrator, even for long term integrations. 

This work shows that if one generalizes the decomposition of e''^"^+^^ to a sum of products, then for a special choice 
of the basis product, the problem reduces to that of simple extrapolation and can be solved easily to any even order. 
When applied to classical Hamiltonian systems, the resulting algorithms are no longer symplectic, but correspond to 
Runge-Kutta-Nystrom (RKN) type integrators. Blanes, Casas and Ros^ and Blanes and Casas^^ have previously 
shown that these extrapolated symplectic integrator enjoys better conservation properties than conventional RKN 
algorithms. Interestingly, multi-product splitting not only reproduces Nystrom integrators at orders p — A and 6 but 
also explain why they only need p —1 force evalutions. At even higher order, these extrapolated RKN algorithms are 
surprisingly competitive with symplectic and other RKN integrators found in the literature. 

II. MULTI-PRODUCT EXPANSION 

If e'^i^+B) were to be decomposed into a sum of products 

^h(A+B) = ^ Cfe [] ea^.^hA^b^,^hB (3) 
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then there is tremendous freedom in the choice of {ck,ak,i,bk,i}- However, for the most general of appHcations, 
including solving time-irreversible diffusion type equationsilii^, one must keep the coefRcients {ak,i,bk i} positive. 
This crucial and deliberate choice distinguishes this work from previous extrapolations, which made no distinction 
between solving time-reversible and time-irreversible problems. If {ak,i,bk^i} were to be positive, then by Sheng's 
theorem^^, any such individual product can at most be second order and some coefficients Ck must be negative. 
Thus without loss of generality, one can choose the simplest basis product as either one of the following symmetric 
second-order splittings 



or 



= e^'^^e'^-^e^''^. 



(4) 
(5) 



The choice of a symmetric product is important, because it has only odd powers of h, 



T2{h) = exp(/i(A + B) + h^Es + h^E^ 



(6) 



where Ei are higher order error commutators of A and B. It then follows that the fcth power of T2 at step size h/k is 
exactly given by 



( ^ j = cxp{h{A + B) + k-^h^E3 + k-'^h^E5 + •••). 



(7) 



Thus the set of powers T^ih/k) has well known error structure and forms a suitable basis for expanding e''("^+^^. For 
example, any two terms with k = I and k — m can approximate e''^"^^^-' to fourth-order via 



c„r2™ - ]+e5ih-'E^) 



(8) 



with obvious solutions 



ci 



m 



m 



(9) 



and error coefficient 



65 = 



(10) 



The RHS of (El) can also be written as 



{l/m) 



1 



(11) 



which coincide with the diagonal elements of the Richardson-Aitken-Neville extrapolation^^ table. Here, we do not 
do the extrapolation numerically, intead, we give an analytical formula for extrapolating to any even order. More 
generally, for a given set of n distinct whole numbers {ki, k2, ---kn}, one can form a 2n-order approximation of e'^i^+B) 
via 



Ji{A+B) 



= ^C,r2^-' (-) +e2„+l(/l'"+^i?2n+l). 



(12) 



The expansion coefficients Ci are determined by a specially simple Vandermonde equation. 



f ^1^ 


1 


- ] 




/Cl\ 

C2 
C3 




= 


,-2(n-l) , -2(n-l) 
'■ I 2 


i,-2("-l) 






\Cn/ 





(13) 
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with closed form solutions 



and error coefRcient, 



n (14) 



e2„+i = (-ir-^n^. (15) 

The closed forms (fH]) and are the key results of this multi-product expansion. All error terms of the same order, 
rather then be forced to zero individually, are extrapolated to zero simultaneously. The Vandermonde equation p3p 
is a special case of extrapolated symplectic algorithms considered by Blanes, Casas and Ros'^*. However, they did not 
obtain the exact solution given by HH). The proof that is the exact solution is given in Appendix A. According to 
Blanes, Casas and Ros^ this class of extrapolated integrator is symplectic to order 2n+3, which at higher orders, can 
easily exceed one's machine precision. That is, at sufficiently high orders, these extrapolated symplectic algorithms 
are, up to machine precision, indistinguishable from truly symplectic algorithms. 

Since a 2n-order expansion is completely characterized by a set of n whole numbers {fci,/c2, •■•fcn}, the resulting 
algorithm will be referred to as the {ki, k2, ...A;„}-integrator. The error coefficient p5|l imphes that, if the number of 
application of 72 is fixed. 



= K (16) 



then the optimal algorithm with the least error is given by a set of n distinct whole numbers {ki} closest to INT (K/n) 
satisfying (|16p . On the other hand, at order 2n, the minimum number of T2 required is given by the natural sequence 
{ki} = {1, 2, 3 ... n}, corresponding to approximating e'^i^^B) ^-^^ 

eHA+B) ^J2^^^j-k +e2„+i(/.2"+i£;2„+i) (17) 

with n(n + l)/2 applications of 72 and an error coefficient of 

e2„+i = (-l)"-^^. (18) 

The simplicity of these results make multi- product splitting algorithms extremely easy to analyze and implement. 
From (jl4|) . the minimum- 72 expansions of orders four to ten are given by 

%ih) = -^-T2ih) + ^Ti(^'^^ (19) 

T(h)- ^ r(h) 64 2 /^M, 6561 //A 16384 390625 

In contrast to a single product decomposition, whose coefficients {ai,bi} are generally irrational and must be deter- 
mined order-by-order numerically with limited precision, the multi-product expansion has only rational coefficients 
and are given analytically by for all even orders. 

The operator extrapolation here is at once more general and simpler than extrapolating solutions of differential 
equations. The general expansion l|17p can be applied to any evolution operator e''''^"'''^-' in terms of 72, with known 
second order error structure ([5]). There is no need to devise and proved, a second-order solution in advance. The low 
order operator extrapolation (I19p has been used previousl y^^'^^ in different contexts. Here, we provide a systematic 
expansion to any even order and point out the connection between symplectic and RKN integrators. 
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III. DERIVING RKN INTEGRATORS 



As emphasized in the last section, the muhi-product expansion (|17[) is an extrapolated operator approximation 
to the evolution operator and can be applied to many different types of equations. Its applications to to 

quantuni^^ and stochstio^i evolution equations have already proven to be highly successful. Here, we will concentrate 
on its use in solving classical dynamic problems. 

If qH^+b) jg decomposed into a sum of products as in ([TT]), then the resulting algorithm will no longer be symplectic. 
This is because for more than one product, cross terms will spoil the usual proof of symplecticity (See Ref.-A, Theorem 
1.). As we will show, multi-product splitting now produces RKN type integrators. 

For solving Hamilton's equation in the form, 



the operators are 



i i 

and where we have abbreviated v — p/m and a(q) = F(q)/m. Thus ^ corresponds to the velocity-form of the 
Verlet algorithm (VV): 

/ ^ 

Vl = vo -I- - a(qo) 
qi = qo + hvi 

V2 = vi + -a(qi) (25) 

and ([5]) corresponds to the position- form of the Verlet algorithm (PV): 

h 

qi = qo + 2 '^o 

Vl = vo /la(qi) 

h , , 

q2 qi + 2 '^i- (26) 

The last numbered variables are the updated variables. Both VV and PV are second order symplectic integrators. 
The operator 7^'^(/i/fc) simply iterating either algorithm fc times at step size hjk and is therefore also symplectic. The 
PV algorithm uses only one force evaluation per update. If 72 is taken to be PV, then the extrapolated integrators 
of order 2n, as examplified by (|19|) - (|22p . only require n(n -t- 1)/2 force-evaluations. The minimal integrator of orders 
4, 6, 8 and 10 requires only 3, 6, 10, and 15 force evaluations, respectively. 

If 72 is taken to the VV algorithm, since all products use the same starting force, each 7^''(/i/fc) for fc > 2 also uses 
only fc force-evaluations. However, 72 (/i) must then evaluate the initial force for the rest of the product and its own 
final force. Thus 72 (/i) requires two force-evaluations, giving the final force count as n(n -|- 1)/2 -|- 1. 

There are other possible extrapolation methods, such as 

T,,«.-lT,(M + HT^(i) (27) 

which is similar to the triplet concatenation for symplectic integrators^i^i^ii^. However, such an iterative extrapolation, 
which triples the number of force-evaluations in going from order 2n to 2n -I- 2, is not competitive with (|17p 's linear 
increase of only n + 1 additional force-evaluations. 

To see how the expansion (IT71) results in RKN integrators, let's takes 7^ to be VV, then the extrapolation (fT5|) 
produces the following fourth-order, {1, 2}-integrator 

q = qo + Vo + Y (ao + 2ai/2) 

v = vo + -^ (ao-F4ai/2 + 2a2/2-ai/i) (28) 
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where we have systematically denoted the force evaluation point of 7^'^(/i/fc) at the intermediate step {i/k)h as q^/j; 
and the resulting force as aj/j, = a(qj/fc), with ap — a(qo). The final force evaluation point of 7^(/i), the midpoint of 
7^'^(/i/2), and the final position of 7^'^(/i/2) are given by 



2 



qi/i = qo + /ivo + yao 



h 

qi/2 = qo + ■^'^o + yao 

q2/2 = qo + /jvo + y (ao + ai/2) . (29) 

The resulting integrator ([28]) appears to require four force-evaluations, in accordance with the formula n{n+ l)/2 + 1. 
However, a key observation here is that force subtractions at the same time step can be combined into a single force 
evaluation. Let 

= q.2/2 - qi/i = —{^1/2 - ao) « 0{h^) (30) 



then the subtraction of two forces gives, 



2a2/2-ai/i 2a(qi/i +^q2) -a(qi/i) (31) 
= a(qi/i + 2<5q2) + 0(/i^) 



1-.2. 



= a ( qo + /ivo + -/i'ai/2 ) + Oih'). (32) 



The total number of force evaluations is reduced from four to three and ([28|) reproduces Nystrom's original fourth- 
order integrato r^^i^^ . Thus we have shown analytically that, the subtraction of two symplectic integrators reproduces 
Nystrom's original integrator. This connection between extrapolated symplectic integrators and Nystrom's integrators 
has not been appreciated previously. 

If T2{h) is taken to be PV, the corresponding {1, 2}-integrator with three force-evaluations is given by 

X = xo -I- /i vo -f y (3ai/4 - ai/2 + a.z/i) 

V = vo + ^ (2ai/4 - ai/2 + 2a3/4) . (33) 

Since PV and VV based algorithms have different force-evaluation points, to avoid confusion, we will denote the 
positions of PV-based integrator as x^/fc. Here, the force-evaluation points are 

h 

Xi/4 = Xo + -Vo 

h 

Xl/2 = Xo + -Vq 

3 h? 

X3/4 = Xo + -/ivo + yai/4. (34) 

All such PV-based integrators are non-FSAL (First Step As Last) RKN integrators in that the force is never evaluated 
initially. The explicit force subtraction will be more susceptible to round-off errors. However, this fourth-order 
integrator is superior to the Nystrom integrator when solving Kepler's orbit and concerns for round-off errors are 
moderated by the wide use of quadruple and multi-precision packages. 

To compare the long time integration error of numerical algorithms, it is useful to solve Kepler's orbit in two 
dimension, 



g3 q3 



with initial conditions 



^ '"y ^ \ " 9x = 1 + e % = (36) 
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where e is the eccentricity of the orbit. The resulting energy and period are respectively —1/2 and 2-k. This set of 
initial conditions is chosen (instead of the usual one in the literature) because it is non-singular as e — > 1 and the 
orbit's semi-major axis is always along the x-axis for all values of e. Both are useful for computing the precession 
error of Kepler's orbit at high eccentricity. 

To gauge the accuracy of the trajectory intrinsically, without any external comparison, we monitor the irreversible 
phase errori^Si^ by computing the rotation angle per period, A^, of the Laplace-Runge-Lenz (LRL) vecto r'^^i'^^ : 

A = V X L - ^ (37) 

with L = q X V. If the orbit is exact, then the LRL vector would remain constant pointing along the semi-major axis 
of the orbit. If A6' 7^ 0, then the orbit will have precessed an angle of A6' after each period. (This is the advance of 
the perihelion of solar planets.) Thus A0 is a much more direct measure of the trajectory's accuracy than the energy 
error. This precession error grows linearly with time as rnA0, where m is the number of periods, even for symplectic 
integrators^Si^i. The corresponding precession error coefficient is extracted by dividing A^ by /i^ using smaller and 
smaller h until 

A0 

lim — = ep (38) 

is independent of h. The coefficient ep as a function of the eccentricity is therefore a unique error signature of each 
fourth-order integrator independent of the step size. 

In FiglU we compare the precession error coefficient of four distinct fourth-order integrators with only three force 
evaluations. The symplectic Forest-Ruth (PR) integrator^, the forward integrator A— 1^, the Nystrom integrator 
(N) given by (|28p . ([5^ and the minimal fourth-ordr {1, 2}-integrator M4, given by ([55]) . The forward integrator A' 
has only positive splitting coefficients and 6^. The convergence of ([55]) is seen near h = 27r/3000; the final value 
used IS h = 27r/5000. The coefficient ep grows steeply with e, by four orders of magnitude from e = 0.4 to e = 0.9, 
but all integrators showed similar dependence on e. It is well known that the precession error of the symplectic FR 
integrator is greater than that of RKN algorithms^Si^ such as N, but it was not known previously that M4 is so much 
better than FR and N. The error coefficient ep at e = 0.9 for FR, N, A' and M4 are -23.1x10^^, 7.1x10'', -1.4x10"* 
and -1.1x10"*, respectively. 



IV. SIXTH ORDER RKN INTEGRATORS 



If T2{h) is taken to be VV, the difference between the two force evaluation points at the same intermediate time 
step, such as (PU)) . is always of order 0{h^) and the resulting force subtration, such as ([5^ . is of order 0{h^). This 
immediately suggests that the three final force evaluations in the {1, 2, 3}-integrator can be collapsed into one, yielding 
a sixth order algorithm with only five force-evaluations. This is indeed the case. The {1, 2, 3}-integrator, according 
to (|20)). is given by 

q = qo + /ivo + Y^(llao + 54ai/3 - 32ai/2 + 27a2/3), (39) 

v = vo + ^(22ao + 162ai/3 - 128ai/2 + 162a2/3 

+ 81a3/3-64a2/2 + 5ai/i), (40) 

with three additional force evaluations at: 

qi/3 = qo + g vo + — ao, (41) 
2 h'^ 

q2/3 = qo + ■::hvQ + — (ao + ai/3) , (42) 
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/j2 

q3/3 = qo + /ivo + — (3ao + 4ai/3 + 2a2/3) . (43) 



Let 



<5q3 = q3/3 - qi/i = l^^i^ai/s + a^^lz ~ 3ao) « 0{h^), (44) 



then the three final-step force-evaluations can be collapsed into one, 

81a3/3 - 64a2/2 + 5ai/i = 22a(qi) + 0{h''), (45) 

with 

qi = qi/i + ^(81^3 - 64(5q2) 

= qo + /ivo + — (l8ai/3 - 16ai/2 + 9a2/3) . (46) 
The force consolidation ([^5]) now renders (1^0]) symmetric, 

V = vo + ^ (22ao + 162ai/3 - 128ai/2 + 162a2/3 + 22Si) . (47) 

This sixth order integrator seems not to be known prior to this work. The well known sixth order integrator with five 
force evaluations due to Albrecht^Si, can be derived from an alternative expansion, 

r«(.,^lr,.,-lr^(i).|^.(i), m 

corresponding to the integrator {1, 2, 4}, 

q = qo + /ivo + ^(7ao + 24ai/4 + 16a2/4 - lOai/2 + 8a3/4), (49) 

h , 

V = vo + — (7ao + 32ai/4 + 32a2/4 - 20ai/2 -I- 32a3/4 

-I- 16a4/4 - 10a2/2 + ai/i), (50) 

with four additional force evaluations at 

h 

qi/4 = qo + ^vo -t- — ao, (51) 
h I? 

(12/4 = qo + -^^0 + Ye + ^1/4) ' (52) 

3 h"^ 

q3/4 = qo + ^^vq + — (3ao + 4ai/4 -t- 2a2/4) , (53) 

q4/4 = qo + /ivo + — (2ao + 3ai/4 + 2a2/4 + a3/4) . (54) 

This algorithm nominally requires eight force evaluations, but forces can now be collapsed at h/2 and at h. At h/2, 
one has 

16a2/4 - lOai/2 = 6a(qt/2) + 0{h<'), (55) 
where the new half time-step evaluation point is 

qt/2 = qo + ^vo + ^ (4ai/4 - ao) . (56) 
Similarly, the three force evaluations at h can be collapsed into one, 

16a4/4 - 10a2/2 + ai/i = 7a(ql) + Oih"), (57) 

where 

ql = qo + h-vo + (6ai/4 -I- 4a2/4 - 5ai/2 + 2a3/4) . (58) 
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The number of force-evaluations can be reduced from eight to five if both q2/4 and qi/2 can be replaced everywhere 
by q*/2- However, replacing a2/4 in q3/4 by a*^2i would yield 

3 /iV 
q3/4 qo + |/ivo + — (^3ao + 4ai/4 + 2a*/2 

( , \ 

- q3/4 + (^^1/2 - ^2/4 j 

= q3/4 + 0(/i5), (59) 

and the resulting 32a*y^ term in (|50p would produce an error of 0{h^), spoiling the algorithm. Remarkably, replacing 
q2/4 and qi/2 everywhere by q^y2 l|58|) defines 

2 



qi = qo + /ivo + ^ (^6ai/4 - a*/2 + 2a5/4 

/j2 

= ^1 " + 4a2/4 - 5ai/2 

= ql + 0(h''), (60) 
and the resulting 0{h^) error of 7a| in ([50| exactly cancels that of 32a3^4! Thus, one recovers Albrecht's integrator: 

q = qo + /ivo + — 7ao + 24ai/4 + 6a*/2 + 8a*/4 +0{h^) (61) 



V - vo + ^ (7ao + 32ai/4 + 12a*/2 + 32a*/4 + 7aj) + 0(/i^). (62) 

Within the context of multi-product expansion, these are the only two sixth-order integrators possible with only 
five force evaluations. Again, we have shown analytically that the extrapolation of three second-order symplectic 
integrators produces sixth-order RKN integrators. 

To compare integrators of the same order but with different number of force evaluations, we compute the precession 
error at e = 0.9 as a function of iV, the number of force-evaluations used. The results for four sixth-order 
algorithms are shown in FigO Y6 is Yoshida's symplectic algorithm^- with the minimum 7 force evaluations. KL6 is 
a much improved version by Kahan and Li^° with 9 force evaluations. A6 and M6 are Albrecht's and the PV-based, 
minimal {l,2,3}-integrator with 5 and 6 force-evaluations respectively. While KL6's error is only half of Y6's at 10^ 
force-evaluations, their respective error are nearly 50 and 100 times larger than those of RKN integrators A6 and M6. 
This calculation was carried out in quadruple-precision using the widely available free software by Miller—. 



V. HIGHER ORDER COMPARISONS 



Since the force reduction process considered is only of 0{h^), this process will no longer be effective for integrators 
beyond sixth-order. Thus beyond sixth-order, some other force reduction processes must be found if the number 
of force-evaluation is to be less than n(n -I- l)/2. At higher orders, since the expansion coefficients are known 
explicitly, it is trivial to write a subroutine for (j26p and simply call it according the multi-product expansion (jl7p . 
Alternatively, one can also construct them from the Aitken-Neville table^^, as done conventionally. It is of interest 
to compare these extrapolated RKN integrators at high orders with existing symplectic and specially derived RKN 
integrators found in the literature. 

In FiglSl we compare some well known higher order integrators with the PV-based, minimal extrapolated integrators 
^T7\ . KL8 and SSIO are eighth and tenth order symplectic integrators by Kahan and Li — and Sofroniou and Spalletta^— 
with 17 and 35 force evaluations respectively. Both are recommended by Ref.?.. Because of their large number of 
force-evaluations, they are not competitive^ in a precision-effort comparison as shown in FigO At iV = 10^, the 
error of KL8 is more than 300 times that of the PV-based {1, 2, 3, 4}-integrator denoted simply as 8 rather than M8. 
Similarly, the error of SSIO at iV = 10^ is about 100 times that of the minimal {1, 2, 3, 4, 5}-integrator denoted as 10. 
DPIO and DP12 are the RKN integrator-pair by Dormand, El-Mikkawy and Prince^ with 17 force evaluations. The 
coefficients for this integrators-pair were taken from Brankin et alM- and converted to quadruple precision. DPIO is 
comparable to 10, which uses 15 force evaluations. DP12's error at TV = 10^ is 10 times smaller than 12. However, 
DP12's superior performance can be quickly matched by going to integrator 14, whose error is then 100 times smaller 
than DP12 at iV = 10^ 
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Fig [3] implies that the most efficient integrator for attaining a given level of precision is never a fixed order integrator. 
As one demands greater and greater precision, one must increase the order of the integrator accordingly. In FigO 
the error curve of the 2n + 2 integrator is plotted only when it is below that of the 2n integrator. The enveloping 
error curve is steeper than any fixed order integrator. This maximum efficiency can be achieved only when one has 
the mean of producing an arbitrary high-order integrator at will, as it is done here. The highest order considered in 
Figl3]is 16 because quadruple precision is inadequate for precision below 10"^''. 



VI. SUMMARY AND CONCLUSIONS 



In this work, we have generalized the single product decomposition of the evolution operator e''^'^"''-^) to that of a 
sum of products. By using powers of T2 as basis products, the multi-product expansion reduces to that of a simple 
extrapolation, with analytically known solutions. Because the extrapolation is formulated on the level of operators, 
it can be applied with great generality in solving diverse evolution equations not only in classical mechanics but also 
in quantum and stochastic dynamics. 

A multi-product expansion loses some desirable properties, such as not being symplectic, unitary, norm preserving, 
etc., as compared to a single product decomposition. However, at high orders, it gains in economy, needing only a 
quadratic, rather than an exponential number of operators. This is a vital consideration when high precision and high 
order results are needed. 

In applying to classical Hamiltonian dynamics, this work showed that extrapolating symplectic integrators produces 
RKN integrators. The process of force consolidation provided a simple explanation for why p—1 force evaluations are 
sufficient for RKN integrators at orders p — A and 6. An order-barrier then naturally arises when forces can longer be 
consolidated at higher p. Moreover, this work showed that, there exists a sequence of easily implemented, 2n-order 
extrapolated RKN integrators, with only n{n + l)/2 force evaluations. Thus any other RKN integrator with fewer 
force evaluations is special and must embody some unique force consolidation schemes. These extrapolated RKN 
integrators are surprisingly competitive with existing symplectic and specially devised RKN algorithms in solving the 
Kepler problem. Since this series of extrapolated RKN integrators can be easily produced by anyone, they can serve 
as a common branchmark by which more sophisticated high order integrators can be judged. For example, of all the 
integrators compared in this work, only DP12 has outperformed its corresponding extrapolated RKN integrator by 
having fewer force evaluations (17) than of n{n + l)/2 (=21). Thus DP12 must possess some highly nontrivial force 
consolidation schemes. 



APPENDIX A: PROVING THE SOLUTION 



To solve the Vandermonde equation for {c,}, let Xi — k^"^ so that (fT5|) reads conventionally 



/ 1 



Xi X2 



1 



(n-l) 



\ 


/Cl\ 






C2 







C3 


= 


/ 


\c„/ 


Vo7 



(Al) 



Consider the usual Lagrange interpolation at n points {xi,X2, • . • , Xn} with values {yi, 2/2, • ■ • , Un}- The interpolating 
n — I degree polynomial is given by 



fix) ^^ViUi^ 

i=l 

where Li(x) are the Lagrange polynomials given by 



(A2) 



Hx)= n 



Since by construction 



{X - Xj) 

Li{xk) = Sik 



(A3) 



(A4) 



10 



the interpolation polynomial (|A2p correctly gives, 

n n 
1=1 i=l 

Now let yi — a;™ for < m < n — 1, then the interpolating polynomial 

n 
i=l 

and the function 

g{x) = x^ 



(A5) 



(A6) 



(A7) 



both interpolate the same set of points. Since interpolating polynomials of the same order are unique, we must have 
f{x) = g{x) and hence 



Y,x^^U{x)^x- 



(A8) 



Writing this out in matrix form yields 



1 1 

Xi X2 
^2 ^2 



X3 



\ ^("-1) ^("-1) 



... 1 X 

Xfi 

x^ 

Xn 




/Li{x)\ 

L2{X) 

Lsix) 


1 ' \ 

X 

= x^ 






\Ln{x)l 





Compare this to (jAip . we immediately see that the solution is 

= L,(0) 



n 



This proof is a simple application of the more general result of inverting the Vandermonde matrijs^. 
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FIG. 1: The orbital precession error coefiicient as a function of the eccentricity e of the Kepler orbit for four fourth-order 
integrators which use only three force evaluations per update. FR is the Forest-Ruth symplectic integrator, N is the original 
Nystrom integrator, A' is a forward integrator and M4 is the multi- product splitting integrator (|33p . 
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FIG. 2: The orbital precession error of sixth-order integrators as a function of of the number of force evaluation A'^. Y6, KL6 
and A6 are integrator of Yoshida, Kahan-Li and Albrecht respectively. M6 is the minimal extrapolated integrator corresponding 

to (unj. 
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FIG. 3: The orbital precession error as a function of the number of force evaluation A^. KL8 and SSIO are Kahan-Li's eighth- 
order and Sofroniou-Spalletta's tenth-order symplectic integrators respectively. DPIO and DP12 are the 10th and 12th order 
RKN pair of Dormand and Prince. The numbers 4, 6, 8, etc., denote the 2n-order integrators corresponding to the minimal 
multi-product expansion (|17p . 



